randn('state',32);
rand('state',32);
load oillogitdata;
id=data(:,1);
s_jt=data(:,2);

x1=[data(:,3:33)];
clear data;
load oilblpiv;
iv=data(:,2:24);
brand=floor(id/100000);
company=floor(brand/10000);

pre_merger_dum=cr_dum(company);
own_dum=pre_merger_dum(1:7,:);


%Base case: 
load oiloutside;
%2.67 q per month
%load oiloutside1;

%load oiloutside2;
%load oiloutside3;
meanprice=data(:,1);
meanshare=data(:,2);
meany=log(reshape(meanshare,7,10))-log(ones(7,10)-kron(ones(7,1),sum(reshape(meanshare,7,10))));


%IV has columns 1 to 9 containing the average price across each region excluding that corresponding to the current observation in each of the 9 time periods.  
%Columns 10 to 15 contain brand dummies
%Column 16 to 24 contain time dummies

%IV=[iv(:,2:2) x1(:,2:14)];
IV=[iv x1(:,4:31)];
clear iv data;
horz=['OLS    IV'];
vert=['constant  ';
      'price     ';
      'private   ';];
%10 REGIONS...7 brands...24 months...
nmkt =24*10;     % number of markets = (# of cities)*(# of months)  %
nbrn =7;     % number of brands per market. if the numebr differs by market this requires some "accounting" vector %
cdid=kron([1:nmkt]',ones(nbrn,1));
cdindex = [nbrn:nbrn:nbrn*nmkt]';       
% create weight matrix
invA = inv([IV'*IV]);



% compute the outside good market share by market
temp = cumsum(s_jt);
sum1 = temp(cdindex,:);
sum1(2:size(sum1,1),:) = diff(sum1);
outshr = 1.0 - sum1(cdid,:);
outshr=outshr.*(outshr>=0);
y = log(s_jt) - log(outshr);

%grand means for elasticity calculations
mean_share=mean(reshape(s_jt,7,24*10)')';
mean_price=mean(reshape(x1(:,1),7,24*10)')';

%ols logit
beta_ols = inv(x1'*x1)*x1'*y;
res = y - x1*beta_ols;
%vcov_ols = 1/size(x1,1)*res'*res*inv(x1'*x1);
vcov_ols = inv(x1'*x1)*(x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1))*inv(x1'*x1);
se_vcov_ols=sqrt(diag(vcov_ols));

%Newey S.E.'s
%Need to get lagged value of xe

x1_lag1=x1(1:1680-70,:);
x1_lag2=x1(1:1680-70*2,:);
x1_lag3=x1(1:1680-70*3,:);


x1_temp1=x1(70+1:1680,:);
x1_temp2=x1(70*2+1:1680,:);
x1_temp3=x1(70*3+1:1680,:);


res_lag1=res(1:1680-70,:);
res_lag2=res(1:1680-70*2,:);
res_lag3=res(1:1680-70*3,:);

res_temp1=res(70+1:1680,:);
res_temp2=res(70*2+1:1680,:);
res_temp3=res(70*3+1:1680,:);


G_0=x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1);
G_1=x1_lag1'.*(ones(size(x1_lag1,2),1)*res_lag1')*(res_temp1*ones(size(x1_temp1,2),1)'.*x1_temp1);
%Need to get twice lagged value of xe
G_2=x1_lag2'.*(ones(size(x1_lag2,2),1)*res_lag2')*(res_temp2*ones(size(x1_temp2,2),1)'.*x1_temp2);
%Need to get three times lagged value of xe
G_3=x1_lag3'.*(ones(size(x1_lag3,2),1)*res_lag3')*(res_temp3*ones(size(x1_temp3,2),1)'.*x1_temp3);



vcov_ols_neweywest = inv(x1'*x1)*(G_0+.5*(G_1+G_1'))*inv(x1'*x1);
se_vcov_ols_newey=sqrt(diag(vcov_ols_neweywest));
clear x1_temp1 x1_temp2 x1_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 G_0 G_1 G_2 G_3;



%ols elasticities at grand means
alpha=beta_ols(1);
o1 = alpha*mean_share*mean_share';
o2 = diag(alpha*mean_share);
omega=(o2-o1);
olselas=omega.*(1./mean_share*mean_price');




% IV regressions
beta_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*IV'*y;
res = y - x1*beta_IV;
IVres = IV.*(res*ones(1,size(IV,2)));
b = IVres'*IVres;
vcov_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*b*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV=sqrt(diag(vcov_IV));


%Newey West S.E.'s here

IV_lag1=IV(1:1680-70,:);
IV_lag2=IV(1:1680-70*2,:);
IV_lag3=IV(1:1680-70*3,:);


IV_temp1=IV(70+1:1680,:);
IV_temp2=IV(70+1:1680,:);
IV_temp3=IV(70+1:1680,:);


res_lag1=res(1:1680-70,:);
res_lag2=res(1:1680-70*2,:);
res_lag3=res(1:1680-70*3,:);

res_temp1=res(70+1:1680,:);
res_temp2=res(70+1:1680,:);
res_temp3=res(70+1:1680,:);

IVres_lag1 = IV.*(res*ones(1,size(IV,2)));
IVres_lag2 = IV.*(res*ones(1,size(IV,2)));
IVres_lag3=IV_lag3.*(res_lag3*ones(1,size(IV_lag3,2)));

G_0=IVres'*IVres;
G_1=IV_lag1'.*(ones(size(IV_lag1,2),1)*res_lag1')*(res_temp1*ones(size(IV_temp1,2),1)'.*IV_temp1);
vcov_IV_neweywest = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*(G_0+.5*(G_1+G_1'))*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV_neweywest=sqrt(diag(vcov_IV_neweywest));

clear IV_lag1 IV_lag2 IV_lag3 IV_temp1 IV_temp2 IV_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 IVres_lag1 IVres_lag2 IVres_lag3;



%iv elasticities at grand means
alpha=beta_IV(1);
o1 = alpha*mean_share*mean_share';
o2 = diag(alpha*mean_share);
omega=(o2-o1);
ivelas=omega.*(1./mean_share*mean_price');



%Approximate first-stage F-stat, leave out vcov d.f. adjustment as N is large relative to number of regresors...
% b_1st=inv(IV'*IV)*IV'*x1(:,1);
% res=x1(:,1)-IV*b_1st;
% vcov_1st=inv(IV'*IV)*(IV'.*(ones(size(IV,2),1)*res')*(res*ones(size(IV,2),1)'.*IV))*inv(IV'*IV);
% se_vcov_1st=sqrt(diag(vcov_1st));
% R=[eye(9) zeros(9,13)];
% q=9;
% F=(R*b_1st)'*inv(R*vcov_1st*R')*(R*b_1st)/q;

alpha=beta_ols(1);
pre_merger_dum=cr_dum(company);
for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
    	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        elas((t-1)*7+1:t*7,:)=omega((t-1)*7+1:t*7,:).*(1./shares*mp');
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==6))+(company==6)*4;
post_dum=cr_dum(company_post);
%post_dum=cr_dum(company);
options=optimset('Tolfun',1e-10,'TolX',1e-10,'Display','on');
first=1;
for i=1:10
	last=i*7;
    expall=exp(meany(:,i)-alpha*meanprice(first:last));
	p_star(first:last,1)=fsolve(@logitbert_eq,meanprice(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end
change=median(reshape(100*((p_star-meanprice)./meanprice),7,10)')

mcflag=zeros(10,1);
first=1;
for i=1:10
    last=i*7;
    expall=exp(meany(:,i)-alpha*meanprice(first:last));
    post_p=meanprice(first:last).*[1+.0141;1-0.2060;1-0.0236;1+0.0396;1-0.0480;1+0.0549;1-0.0448];
    [mc_post(first:last,1),y,mcflag(i,1)]=fsolve(@mc_logitbert,mc_hat(first:last),options,expall,alpha,post_p,post_dum(first:last,:));
    first=last+1;
end
olsmc_change=reshape(100*((mc_post-mc_hat)./mc_hat),7,10)';
olsmcchange=median(olsmc_change);


alpha=beta_IV(1);
pre_merger_dum=cr_dum(company);
for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
    	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        elas((t-1)*7+1:t*7,:)=omega((t-1)*7+1:t*7,:).*(1./shares*mp');
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==6))+(company==6)*4;
post_dum=cr_dum(company_post);
%post_dum=cr_dum(company);
options=optimset('Tolfun',1e-15,'TolX',1e-15,'Display','on');
first=1;
for i=1:10
	last=i*7;
    expall=exp(meany(:,i)-alpha*meanprice(first:last));
	p_star(first:last,1)=fsolve(@logitbert_eq,meanprice(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end
ivchange=median(reshape(100*((p_star-meanprice)./meanprice),7,10)')

mcflag=zeros(10,1);
first=1;
for i=1:10
    last=i*7;
    expall=exp(meany(:,i)-alpha*meanprice(first:last));
    post_p=meanprice(first:last).*[1+.0141;1-0.2060;1-0.0236;1+0.0396;1-0.0480;1+0.0549;1-0.0448];
    [mc_post(first:last,1),y,mcflag(i,1)]=fsolve(@mc_logitbert,mc_hat(first:last),options,expall,alpha,post_p,post_dum(first:last,:));
    first=last+1;
end
ivmc_change=reshape(100*((mc_post-mc_hat)./mc_hat),7,10)';
ivmcchange=median(ivmc_change);



%number of bootstrap draws
sims=1000;
pchange=zeros(sims,7);
melas=zeros(7,7,sims);

for i=1:sims
    bs_beta_IV=beta_IV+chol(vcov_IV_neweywest)'*randn(size(beta_IV,1),1);
    bs_beta_ols=beta_ols+chol(vcov_ols_neweywest)'*randn(size(beta_ols,1),1);
    alpha=bs_beta_IV(1);         

    markup=zeros(70,1);
    p_star=zeros(70,1);
    % compute demand elasticities at overall means
    o1=alpha*mean_share*mean_share';
    o2=diag(alpha*mean_share);
    omega=(o2-o1);
    melas(:,:,i)=omega.*(1./mean_share*mean_price');
    
for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
     	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==6))+(company==6)*4;
post_dum=cr_dum(company_post);
options=optimset('Tolfun',1e-10,'TolX',1e-14,'Display','on');
first=1;
for j=1:10
	last=j*7;
    expall=exp(meany(:,j)-alpha*meanprice(first:last));
	p_star(first:last,1)=fsolve(@logitbert_eq,meanprice(first:last),options,expall,alpha,mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end

pchange(i,:)=median(reshape(100*((p_star-meanprice)./meanprice),7,10)');
   

end

